{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 3     \n",
    "\n",
    "Use Tablesaw to plot graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load ../../utils/djl-imports\n",
    "%load ../../utils/plot-utils\n",
    "%load ../../utils/Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import org.apache.commons.lang3.ArrayUtils;\n",
    "\n",
    "NDManager manager = NDManager.newBaseManager();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Logistic Regression for Binary Classification\n",
    "\n",
    "In multiclass classification we typically use the exponential model \n",
    "\n",
    "$$p(y|\\mathbf{o}) = \\mathrm{softmax}(\\mathbf{o})_y = \\frac{\\exp(o_y)}{\\sum_{y'} \\exp(o_{y'})}$$\n",
    "\n",
    "1.1. Show that this parametrization has a spurious degree of freedom. That is, show that both $\\mathbf{o}$ and $\\mathbf{o} + c$ with $c \\in \\mathbb{R}$ lead to the same probability estimate.\n",
    "1.2. For binary classification, i.e. whenever we have only two classes $\\{-1, 1\\}$, we can arbitrarily set $o_{-1} = 0$. Using the shorthand $o = o_1$ show that this is equivalent to \n",
    "\n",
    "$$p(y=1|o) = \\frac{1}{1 + \\exp(-o)}$$\n",
    "\n",
    "1.3. Show that the log-likelihood loss (often called logistic loss) for labels $y \\in \\{-1, 1\\}$ is thus given by \n",
    "\n",
    "$$-\\log p(y|o) = \\log (1 + \\exp(-y \\cdot o))$$\n",
    "\n",
    "1.4. Show that for $y = 1$ the logistic loss asymptotes to $o$ for $o \\to \\infty$ and to $\\exp(o)$ for $o \\to -\\infty$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Logistic Regression and Autograd\n",
    "\n",
    "1. Implement the binary logistic loss $l(y,o) = \\log (1 + \\exp(-y \\cdot o))$ in DJL\n",
    "1. Plot its values for $y \\in \\{-1, 1\\}$ over the range of $o \\in [-5, 5]$.\n",
    "1. Plot its derivative with respect to $o$ for $o \\in [-5, 5]$ using 'GradientCollector'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-07T00:58:19.464756Z",
     "start_time": "2019-02-07T00:58:19.459943Z"
    }
   },
   "outputs": [],
   "source": [
    "float loss(float y, float o) {\n",
    "    // Add your loss function here\n",
    "    float l = 0;\n",
    "    return l;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Ohm's Law\n",
    "\n",
    "Imagine that you're a young physicist, maybe named [Georg Simon Ohm](https://en.wikipedia.org/wiki/Georg_Ohm), trying to figure out how current and voltage depend on each other for resistors. You have some idea but you aren't quite sure yet whether the dependence is linear or quadratic. So you take some measurements, conveniently given to you as 'NDArrays' in DJL. They are indicated by 'current' and 'voltage'.\n",
    "\n",
    "Your goal is to use least mean squares regression to identify the coefficients for the following three models using automatic differentiation and least mean squares regression. The three models are:\n",
    "\n",
    "1. Quadratic model where $\\mathrm{voltage} = c + r \\cdot \\mathrm{current} + q \\cdot \\mathrm{current}^2$.\n",
    "1. Linear model where $\\mathrm{voltage} = c + r \\cdot \\mathrm{current}$.\n",
    "1. Ohm's law where $\\mathrm{voltage} = r \\cdot \\mathrm{current}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-07T00:58:19.477109Z",
     "start_time": "2019-02-07T00:58:19.468122Z"
    }
   },
   "outputs": [],
   "source": [
    "var current = manager.create(new float[]{1.5420291f, 1.8935232f, 2.1603365f, 2.5381863f, 2.893443f,\n",
    "                    3.838855f, 3.92542f, 4.223369f, 4.235571f, 4.273397f,\n",
    "                    4.9332876f, 6.4704757f, 6.517571f, 6.87826f, 7.0009003f, \n",
    "                    7.035741f, 7.278681f, 7.7561755f, 9.121138f, 9.728281f});\n",
    "var voltage = manager.create(new float[]{63.802246f, 80.036026f, 91.4903f, 108.28776f, 122.781975f,\n",
    "                    161.36314f, 166.50816f, 176.16772f, 180.29395f, 179.09758f,\n",
    "                    206.21027f, 272.71857f, 272.24033f, 289.54745f, 293.8488f,\n",
    "                    295.2281f, 306.62274f, 327.93243f, 383.16296f, 408.65967f})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Entropy\n",
    "\n",
    "Let's compute the *binary* entropy of a number of interesting data sources. \n",
    "\n",
    "1. Assume that you're watching the output generated by a [monkey at a typewriter](https://en.wikipedia.org/wiki/File:Chimpanzee_seated_at_typewriter.jpg). The monkey presses any of the $44$ keys of the typewriter at random (you can assume that it has not discovered any special keys or the shift key yet). How many bits of randomness per character do you observe?\n",
    "1. Unhappy with the monkey you replaced it by a drunk typesetter. It is able to generate words, albeit not coherently. Instead, it picks a random word out of a vocabulary of $2,000$ words. Moreover, assume that the average length of a word is $4.5$ letters in English. How many bits of randomness do you observe now?\n",
    "1. Still unhappy with the result you replace the typesetter by a high quality language model. These can obtain perplexity numbers as low as 20 points per character. The perplexity is defined as a length normalized probability, i.e.\n",
    "\n",
    "$$\\mathrm{PPL}(x) = \\left[p(x)\\right]^{1/\\mathrm{length}(x)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Wien's Approximation for the Temperature (bonus)\n",
    "\n",
    "We will now abuse DJL to estimate the temperature of a black body. The energy emanated from a black body is given by Wien's approximation.\n",
    "\n",
    "$$B_\\lambda(T) = \\frac{2 h c^2}{\\lambda^5} \\exp\\left(-\\frac{h c}{\\lambda k T}\\right)$$\n",
    "\n",
    "That is, the amount of energy depends on the fifth power of the wavelength $\\lambda$ and the temperature $T$ of the body. The latter ensures a cutoff beyond a temperature-characteristic peak. Let us define this and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-07T00:58:19.805358Z",
     "start_time": "2019-02-07T00:58:19.479818Z"
    }
   },
   "outputs": [],
   "source": [
    "// Lightspeed\n",
    "float c = 299792458;\n",
    "// Planck's constant\n",
    "float h = (float) 6.62607004e-34;\n",
    "// Boltzmann constant\n",
    "float k = (float) 1.38064852e-23;\n",
    "// Wavelength scale (nanometers)\n",
    "float lamscale = (float) 1e-6;\n",
    "// Pulling out all powers of 10 upfront\n",
    "float pOut = 2 * h * (float) Math.pow(c, 2) / (float) Math.pow(lamscale, 5);\n",
    "float pIn = (h / k) * (c / lamscale);\n",
    "\n",
    "// Wien's law\n",
    "NDArray wien(NDArray lam, float t) {\n",
    "    return lam.pow(-5).mul(pOut).mul(lam.mul(t).pow(-1).mul(-pIn).exp());\n",
    "}\n",
    "\n",
    "// Plot the radiance for a few different temperatures\n",
    "var lam = manager.arange(0, 100, 0.01f);\n",
    "float[] lamArray = lam.toFloatArray();\n",
    "\n",
    "// To hold all data\n",
    "float[] lambdas = new float[0];\n",
    "float[] radiances = new float[0];\n",
    "String[] T = new String[0];\n",
    "\n",
    "for (float t : new float[]{10, 100, 150, 200, 250, 300, 350}) {\n",
    "    float[] radianceArray = wien(lam, t, pOut, pIn).toFloatArray();\n",
    "    lambdas = ArrayUtils.addAll(lambdas, lamArray);\n",
    "    radiances = ArrayUtils.addAll(radiances, radianceArray);\n",
    "    String tString = String.format(\"T=%dK\", (int) t);\n",
    "    String[] tArray = new String[lamArray.length];\n",
    "    Arrays.fill(tArray, tString);\n",
    "    T = ArrayUtils.addAll(T, tArray);\n",
    "}\n",
    "\n",
    "Table data = Table.create(\"data\")\n",
    "    .addColumns(\n",
    "        DoubleColumn.create(\"lambda\", Functions.floatToDoubleArray(lambdas)),\n",
    "        DoubleColumn.create(\"radiance\", Functions.floatToDoubleArray(radiances)),\n",
    "        StringColumn.create(\"T\", T)\n",
    ");\n",
    "LinePlot.create(\"radiance vs. lambda\", data, \"lambda\", \"radiance\", \"T\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we assume that we are a fearless physicist measuring some data. Of course, we need to pretend that we don't really know the temperature. But we measure the radiation at a few wavelengths.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-07T00:58:20.063122Z",
     "start_time": "2019-02-07T00:58:19.808000Z"
    }
   },
   "outputs": [],
   "source": [
    "// Real temperature is approximately 0C\n",
    "float realtemp = 273;\n",
    "// We observe at 3000nm up to 20,000nm wavelength\n",
    "var wavelengths = manager.arange(3,20,2f);\n",
    "float[] wavelengthArray = wavelengths.toFloatArray();\n",
    "// Our infrared filters are pretty lousy ...\n",
    "var delta = manager.randomNormal(new Shape(wavelengths.size()));\n",
    "\n",
    "var radiance = wien(wavelengths.add(delta), realtemp);\n",
    "var radianceTrue = wien(wavelengths, realtemp);\n",
    "\n",
    "float[] wavelengthData = new float[0];\n",
    "float[] radianceData = new float[0];\n",
    "String[] typeData = new String[0];\n",
    "\n",
    "wavelengthData = ArrayUtils.addAll(wavelengthArray, wavelengthArray);\n",
    "radianceData = ArrayUtils.addAll(radiance.toFloatArray(), radianceTrue.toFloatArray());\n",
    "String[] measuredType = new String[wavelengthArray.length];\n",
    "Arrays.fill(measuredType, \"measured\");\n",
    "String[] trueType = new String[wavelengthArray.length];\n",
    "Arrays.fill(trueType, \"true\");\n",
    "typeData = ArrayUtils.addAll(measuredType, trueType);\n",
    "\n",
    "Table data = Table.create(\"data\")\n",
    "    .addColumns(\n",
    "        DoubleColumn.create(\"wavelength\", Functions.floatToDoubleArray(wavelengthData)),\n",
    "        DoubleColumn.create(\"radiance\", Functions.floatToDoubleArray(radianceData)),\n",
    "        StringColumn.create(\"type\", typeData)\n",
    ");\n",
    "LinePlot.create(\"radiance vs. wavelength\", data, \"wavelength\", \"radiance\", \"type\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T19:41:18.204791Z",
     "start_time": "2019-01-29T19:41:18.201091Z"
    }
   },
   "source": [
    "Use DJL to estimate the real temperature based on the variables `wavelengths` and `radiance`. \n",
    "\n",
    "* You can use Wien's law implementation `wien(lam,t)` as your forward model. \n",
    "* Use the loss function $l(y,y') = (\\log y - \\log y')^2$ to measure accuracy."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Java",
   "language": "java",
   "name": "java"
  },
  "language_info": {
   "codemirror_mode": "java",
   "file_extension": ".jshell",
   "mimetype": "text/x-java-source",
   "name": "Java",
   "pygments_lexer": "java",
   "version": "14.0.2+12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
